home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / gnu / calc202a.lha / calc-2.02a / calc-mat.el < prev    next >
Lisp/Scheme  |  1993-06-01  |  10KB  |  379 lines
  1. ;; Calculator for GNU Emacs, part II [calc-mat.el]
  2. ;; Copyright (C) 1990, 1991, 1992, 1993 Free Software Foundation, Inc.
  3. ;; Written by Dave Gillespie, daveg@synaptics.com.
  4.  
  5. ;; This file is part of GNU Emacs.
  6.  
  7. ;; GNU Emacs is distributed in the hope that it will be useful,
  8. ;; but WITHOUT ANY WARRANTY.  No author or distributor
  9. ;; accepts responsibility to anyone for the consequences of using it
  10. ;; or for whether it serves any particular purpose or works at all,
  11. ;; unless he says so in writing.  Refer to the GNU Emacs General Public
  12. ;; License for full details.
  13.  
  14. ;; Everyone is granted permission to copy, modify and redistribute
  15. ;; GNU Emacs, but only under the conditions described in the
  16. ;; GNU Emacs General Public License.   A copy of this license is
  17. ;; supposed to have been given to you along with GNU Emacs so you
  18. ;; can know your rights and responsibilities.  It should be in a
  19. ;; file named COPYING.  Among other things, the copyright notice
  20. ;; and this notice must be preserved on all copies.
  21.  
  22.  
  23.  
  24. ;; This file is autoloaded from calc-ext.el.
  25. (require 'calc-ext)
  26.  
  27. (require 'calc-macs)
  28.  
  29. (defun calc-Need-calc-mat () nil)
  30.  
  31.  
  32. (defun calc-mdet (arg)
  33.   (interactive "P")
  34.   (calc-slow-wrapper
  35.    (calc-unary-op "mdet" 'calcFunc-det arg))
  36. )
  37.  
  38. (defun calc-mtrace (arg)
  39.   (interactive "P")
  40.   (calc-slow-wrapper
  41.    (calc-unary-op "mtr" 'calcFunc-tr arg))
  42. )
  43.  
  44. (defun calc-mlud (arg)
  45.   (interactive "P")
  46.   (calc-slow-wrapper
  47.    (calc-unary-op "mlud" 'calcFunc-lud arg))
  48. )
  49.  
  50.  
  51. ;;; Coerce row vector A to be a matrix.  [V V]
  52. (defun math-row-matrix (a)
  53.   (if (and (Math-vectorp a)
  54.        (not (math-matrixp a)))
  55.       (list 'vec a)
  56.     a)
  57. )
  58.  
  59. ;;; Coerce column vector A to be a matrix.  [V V]
  60. (defun math-col-matrix (a)
  61.   (if (and (Math-vectorp a)
  62.        (not (math-matrixp a)))
  63.       (cons 'vec (mapcar (function (lambda (x) (list 'vec x))) (cdr a)))
  64.     a)
  65. )
  66.  
  67.  
  68.  
  69. ;;; Multiply matrices A and B.  [V V V]
  70. (defun math-mul-mats (a b)
  71.   (let ((mat nil)
  72.     (cols (length (nth 1 b)))
  73.     row col ap bp accum)
  74.     (while (setq a (cdr a))
  75.       (setq col cols
  76.         row nil)
  77.       (while (> (setq col (1- col)) 0)
  78.     (setq ap (cdr (car a))
  79.           bp (cdr b)
  80.           accum (math-mul (car ap) (nth col (car bp))))
  81.     (while (setq ap (cdr ap) bp (cdr bp))
  82.       (setq accum (math-add accum (math-mul (car ap) (nth col (car bp))))))
  83.     (setq row (cons accum row)))
  84.       (setq mat (cons (cons 'vec row) mat)))
  85.     (cons 'vec (nreverse mat)))
  86. )
  87.  
  88. (defun math-mul-mat-vec (a b)
  89.   (cons 'vec (mapcar (function (lambda (row)
  90.                  (math-dot-product row b)))
  91.              (cdr a)))
  92. )
  93.  
  94.  
  95.  
  96. (defun calcFunc-tr (mat)   ; [Public]
  97.   (if (math-square-matrixp mat)
  98.       (math-matrix-trace-step 2 (1- (length mat)) mat (nth 1 (nth 1 mat)))
  99.     (math-reject-arg mat 'square-matrixp))
  100. )
  101.  
  102. (defun math-matrix-trace-step (n size mat sum)
  103.   (if (<= n size)
  104.       (math-matrix-trace-step (1+ n) size mat
  105.                   (math-add sum (nth n (nth n mat))))
  106.     sum)
  107. )
  108.  
  109.  
  110. ;;; Matrix inverse and determinant.
  111. (defun math-matrix-inv-raw (m)
  112.   (let ((n (1- (length m))))
  113.     (if (<= n 3)
  114.     (let ((det (math-det-raw m)))
  115.       (and (not (math-zerop det))
  116.            (math-div
  117.         (cond ((= n 1) 1)
  118.               ((= n 2)
  119.                (list 'vec
  120.                  (list 'vec
  121.                    (nth 2 (nth 2 m))
  122.                    (math-neg (nth 2 (nth 1 m))))
  123.                  (list 'vec
  124.                    (math-neg (nth 1 (nth 2 m)))
  125.                    (nth 1 (nth 1 m)))))
  126.               ((= n 3)
  127.                (list 'vec
  128.                  (list 'vec
  129.                    (math-sub (math-mul (nth 3 (nth 3 m))
  130.                                (nth 2 (nth 2 m)))
  131.                          (math-mul (nth 3 (nth 2 m))
  132.                                (nth 2 (nth 3 m))))
  133.                    (math-sub (math-mul (nth 3 (nth 1 m))
  134.                                (nth 2 (nth 3 m)))
  135.                          (math-mul (nth 3 (nth 3 m))
  136.                                (nth 2 (nth 1 m))))
  137.                    (math-sub (math-mul (nth 3 (nth 2 m))
  138.                                (nth 2 (nth 1 m)))
  139.                          (math-mul (nth 3 (nth 1 m))
  140.                                (nth 2 (nth 2 m)))))
  141.                  (list 'vec
  142.                    (math-sub (math-mul (nth 3 (nth 2 m))
  143.                                (nth 1 (nth 3 m)))
  144.                          (math-mul (nth 3 (nth 3 m))
  145.                                (nth 1 (nth 2 m))))
  146.                    (math-sub (math-mul (nth 3 (nth 3 m))
  147.                                (nth 1 (nth 1 m)))
  148.                          (math-mul (nth 3 (nth 1 m))
  149.                                (nth 1 (nth 3 m))))
  150.                    (math-sub (math-mul (nth 3 (nth 1 m))
  151.                                (nth 1 (nth 2 m)))
  152.                          (math-mul (nth 3 (nth 2 m))
  153.                                (nth 1 (nth 1 m)))))
  154.                  (list 'vec
  155.                    (math-sub (math-mul (nth 2 (nth 3 m))
  156.                                (nth 1 (nth 2 m)))
  157.                          (math-mul (nth 2 (nth 2 m))
  158.                                (nth 1 (nth 3 m))))
  159.                    (math-sub (math-mul (nth 2 (nth 1 m))
  160.                                (nth 1 (nth 3 m)))
  161.                          (math-mul (nth 2 (nth 3 m))
  162.                                (nth 1 (nth 1 m))))
  163.                    (math-sub (math-mul (nth 2 (nth 2 m))
  164.                                (nth 1 (nth 1 m)))
  165.                          (math-mul (nth 2 (nth 1 m))
  166.                                (nth 1 (nth 2 m))))))))
  167.         det)))
  168.       (let ((lud (math-matrix-lud m)))
  169.     (and lud
  170.          (math-lud-solve lud (calcFunc-idn 1 n))))))
  171. )
  172.  
  173. (defun calcFunc-det (m)
  174.   (if (math-square-matrixp m)
  175.       (math-with-extra-prec 2 (math-det-raw m))
  176.     (if (and (eq (car-safe m) 'calcFunc-idn)
  177.          (or (math-zerop (nth 1 m))
  178.          (math-equal-int (nth 1 m) 1)))
  179.     (nth 1 m)
  180.       (math-reject-arg m 'square-matrixp)))
  181. )
  182.  
  183. (defun math-det-raw (m)
  184.   (let ((n (1- (length m))))
  185.     (cond ((= n 1)
  186.        (nth 1 (nth 1 m)))
  187.       ((= n 2)
  188.        (math-sub (math-mul (nth 1 (nth 1 m))
  189.                    (nth 2 (nth 2 m)))
  190.              (math-mul (nth 2 (nth 1 m))
  191.                    (nth 1 (nth 2 m)))))
  192.       ((= n 3)
  193.        (math-sub
  194.         (math-sub
  195.          (math-sub
  196.           (math-add
  197.            (math-add
  198.         (math-mul (nth 1 (nth 1 m))
  199.               (math-mul (nth 2 (nth 2 m))
  200.                     (nth 3 (nth 3 m))))
  201.         (math-mul (nth 2 (nth 1 m))
  202.               (math-mul (nth 3 (nth 2 m))
  203.                     (nth 1 (nth 3 m)))))
  204.            (math-mul (nth 3 (nth 1 m))
  205.              (math-mul (nth 1 (nth 2 m))
  206.                    (nth 2 (nth 3 m)))))
  207.           (math-mul (nth 3 (nth 1 m))
  208.             (math-mul (nth 2 (nth 2 m))
  209.                   (nth 1 (nth 3 m)))))
  210.          (math-mul (nth 1 (nth 1 m))
  211.                (math-mul (nth 3 (nth 2 m))
  212.                  (nth 2 (nth 3 m)))))
  213.         (math-mul (nth 2 (nth 1 m))
  214.               (math-mul (nth 1 (nth 2 m))
  215.                 (nth 3 (nth 3 m))))))
  216.       (t (let ((lud (math-matrix-lud m)))
  217.            (if lud
  218.            (let ((lu (car lud)))
  219.              (math-det-step n (nth 2 lud)))
  220.          0)))))
  221. )
  222.  
  223. (defun math-det-step (n prod)
  224.   (if (> n 0)
  225.       (math-det-step (1- n) (math-mul prod (nth n (nth n lu))))
  226.     prod)
  227. )
  228.  
  229. ;;; This returns a list (LU index d), or NIL if not possible.
  230. ;;; Argument M must be a square matrix.
  231. (defun math-matrix-lud (m)
  232.   (let ((old (assoc m math-lud-cache))
  233.     (context (list calc-internal-prec calc-prefer-frac)))
  234.     (if (and old (equal (nth 1 old) context))
  235.     (cdr (cdr old))
  236.       (let* ((lud (catch 'singular (math-do-matrix-lud m)))
  237.          (entry (cons context lud)))
  238.     (if old
  239.         (setcdr old entry)
  240.       (setq math-lud-cache (cons (cons m entry) math-lud-cache)))
  241.     lud)))
  242. )
  243. (defvar math-lud-cache nil)
  244.  
  245. ;;; Numerical Recipes section 2.3; implicit pivoting omitted.
  246. (defun math-do-matrix-lud (m)
  247.   (let* ((lu (math-copy-matrix m))
  248.      (n (1- (length lu)))
  249.      i (j 1) k imax sum big
  250.      (d 1) (index nil))
  251.     (while (<= j n)
  252.       (setq i 1
  253.         big 0
  254.         imax j)
  255.       (while (< i j)
  256.     (math-working "LUD step" (format "%d/%d" j i))
  257.     (setq sum (nth j (nth i lu))
  258.           k 1)
  259.     (while (< k i)
  260.       (setq sum (math-sub sum (math-mul (nth k (nth i lu))
  261.                         (nth j (nth k lu))))
  262.         k (1+ k)))
  263.     (setcar (nthcdr j (nth i lu)) sum)
  264.     (setq i (1+ i)))
  265.       (while (<= i n)
  266.     (math-working "LUD step" (format "%d/%d" j i))
  267.     (setq sum (nth j (nth i lu))
  268.           k 1)
  269.     (while (< k j)
  270.       (setq sum (math-sub sum (math-mul (nth k (nth i lu))
  271.                         (nth j (nth k lu))))
  272.         k (1+ k)))
  273.     (setcar (nthcdr j (nth i lu)) sum)
  274.     (let ((dum (math-abs-approx sum)))
  275.       (if (Math-lessp big dum)
  276.           (setq big dum
  277.             imax i)))
  278.     (setq i (1+ i)))
  279.       (if (> imax j)
  280.       (setq lu (math-swap-rows lu j imax)
  281.         d (- d)))
  282.       (setq index (cons imax index))
  283.       (let ((pivot (nth j (nth j lu))))
  284.     (if (math-zerop pivot)
  285.         (throw 'singular nil)
  286.       (setq i j)
  287.       (while (<= (setq i (1+ i)) n)
  288.         (setcar (nthcdr j (nth i lu))
  289.             (math-div (nth j (nth i lu)) pivot)))))
  290.       (setq j (1+ j)))
  291.     (list lu (nreverse index) d))
  292. )
  293.  
  294. (defun math-swap-rows (m r1 r2)
  295.   (or (= r1 r2)
  296.       (let* ((r1prev (nthcdr (1- r1) m))
  297.          (row1 (cdr r1prev))
  298.          (r2prev (nthcdr (1- r2) m))
  299.          (row2 (cdr r2prev))
  300.          (r2next (cdr row2)))
  301.     (setcdr r2prev row1)
  302.     (setcdr r1prev row2)
  303.     (setcdr row2 (cdr row1))
  304.     (setcdr row1 r2next)))
  305.   m
  306. )
  307.  
  308.  
  309. (defun math-lud-solve (lud b &optional need)
  310.   (if lud
  311.       (let* ((x (math-copy-matrix b))
  312.          (n (1- (length x)))
  313.          (m (1- (length (nth 1 x))))
  314.          (lu (car lud))
  315.          (col 1)
  316.          i j ip ii index sum)
  317.     (while (<= col m)
  318.       (math-working "LUD solver step" col)
  319.       (setq i 1
  320.         ii nil
  321.         index (nth 1 lud))
  322.       (while (<= i n)
  323.         (setq ip (car index)
  324.           index (cdr index)
  325.           sum (nth col (nth ip x)))
  326.         (setcar (nthcdr col (nth ip x)) (nth col (nth i x)))
  327.         (if (null ii)
  328.         (or (math-zerop sum)
  329.             (setq ii i))
  330.           (setq j ii)
  331.           (while (< j i)
  332.         (setq sum (math-sub sum (math-mul (nth j (nth i lu))
  333.                           (nth col (nth j x))))
  334.               j (1+ j))))
  335.         (setcar (nthcdr col (nth i x)) sum)
  336.         (setq i (1+ i)))
  337.       (while (>= (setq i (1- i)) 1)
  338.         (setq sum (nth col (nth i x))
  339.           j i)
  340.         (while (<= (setq j (1+ j)) n)
  341.           (setq sum (math-sub sum (math-mul (nth j (nth i lu))
  342.                         (nth col (nth j x))))))
  343.         (setcar (nthcdr col (nth i x))
  344.             (math-div sum (nth i (nth i lu)))))
  345.       (setq col (1+ col)))
  346.     x)
  347.     (and need
  348.      (math-reject-arg need "*Singular matrix")))
  349. )
  350.  
  351. (defun calcFunc-lud (m)
  352.   (if (math-square-matrixp m)
  353.       (or (math-with-extra-prec 2
  354.         (let ((lud (math-matrix-lud m)))
  355.           (and lud
  356.            (let* ((lmat (math-copy-matrix (car lud)))
  357.               (umat (math-copy-matrix (car lud)))
  358.               (n (1- (length (car lud))))
  359.               (perm (calcFunc-idn 1 n))
  360.               i (j 1))
  361.              (while (<= j n)
  362.                (setq i 1)
  363.                (while (< i j)
  364.              (setcar (nthcdr j (nth i lmat)) 0)
  365.              (setq i (1+ i)))
  366.                (setcar (nthcdr j (nth j lmat)) 1)
  367.                (while (<= (setq i (1+ i)) n)
  368.              (setcar (nthcdr j (nth i umat)) 0))
  369.                (setq j (1+ j)))
  370.              (while (>= (setq j (1- j)) 1)
  371.                (let ((pos (nth (1- j) (nth 1 lud))))
  372.              (or (= pos j)
  373.                  (setq perm (math-swap-rows perm j pos)))))
  374.              (list 'vec perm lmat umat)))))
  375.       (math-reject-arg m "*Singular matrix"))
  376.     (math-reject-arg m 'square-matrixp))
  377. )
  378.  
  379.